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The structure of a simple liquid may be characterised in terms of ground state clusters of small 
numbers of atoms of that same liquid. Here we use this sensitive structural probe to consider the 
effect of a liquid-vapour interface upon the liquid structure. At higher temperatures (above around 
half the critical temperature) we find that the predominant effect of the interface is to reduce the local 
density, which significantly suppresses the local cluster populations. At lower temperatures, however, 
pronounced interfacial layering is found. This appears to be connected with significant orient at ional 
ordering of clusters based on 3- and 5-membered rings, with the rings aligning perpendicular and 
parallel to the interface respectively. At all temperatures, we find that the population of five-fold 
symmetric structures is suppressed, rather than enhanced, close to the interface. 



I. INTRODUCTION 

What is the structure of the transition zone between 
liquid and vapour? In 1893 van der Waals^ ^ described 
the liquid-vapour interface as a smooth density profile 
p{z)^ shortly preceded by Rayleigh^. 

Early computational^ and theoretical work based on 
the Born-Green- Yvon equation ^ suggested a significant 
degree of surface layering in the Lennard- Jones model of 
simple liquids. However, further simulation^ and den- 
sity functional theory ^ ^ showed that little surface lay- 
ering is expected for materials based on van der Waals 
interactions, and that the interfacial profile is dominated 
by a smooth change in density. Relative to the critical 
temperature Tc, metals typically remain liquid at rather 
lower temperatures than Lennard- Jones like materials, 
and here pronounced surface layering is indeed found^ 
An important related problem to the free interface is the 
confinement introduced by a hard wall, crystal or similar 
external field, which may be tackled with den sity f unc- 
tional theor};^^^, simulation^ and experiment! ^ I ^^ l In 
this case, layering is more generally observed. 

In experiments on conventional materials, the mea- 
surement of (intrinsic) interfacial profiles is complicated 
by collective surface excitations in the form of capillary 
waveJI^HI^. For those systems where surface layering is 
pronounced such as metals, the lengthscale above which 
capillary wave effects become dominant is around 100 
nm. This lengthscale is larger than that typically encoun- 
tered in computational studies, where capillary waves 
only cause a small broadening of surface layering^, but 
cannot be neglected in experiment. Nonetheless surface 
layering has been demonstrated in experiments, for ex- 
ample in liquid M.eicmy^'^ . 

A further class of materials which may be directly in- 
vestigated in real space due to their mesoscopic length- 



scale are colloidal suspensions. Most notably, colloid- 
polymer mixtures exhibit (colloidal) vapour-liquid phase 
separation^^ moreover the hard-sphere interactions 
that closely approximate such systems are amenable to 
theor}/^^-. Here capillary waves have been directly vi- 
sualised^. So far, to our knowledge, interfaces in col- 
loidal systems with long-ranged interactions^^ which can 
be deeply cooled and might exhibit pronounced layering 
(the equivalent of metals) have not been studied. How- 
ever, single-particle resolution enables other approaches 
to be used. Notably, a method introduced by Chacon and 
TarazonsP^ which pins a plane to a layer of surface parti- 
cles indeed revealed density oscillations perpendicular to 
the pinned surface in a colloidal liquid^^. 

Here we consider the perturbation induced by the in- 
terface on the liquid structure at the molecular /atomic 
level. We consider two different liquid systems. The first 
is a truncated and shifted Lennard- Jones (LJ) liquicP^^ 
whose triple temperature T^-/ lies at T^/ /T^'^ = 0.630 
where T^^^ is the critical temperature. To observe a 
stronger effect of the liquid-gas interface, we also consider 
a model potential which approximates Sodium (Na), 
whose triple point lies at a rather lowe r tem perature rel- 
ative to criticality of T^^^/T^^^ = 0.22^ . This model 
is known to exhibit strong surface layering around the 
triple point, which is also seen in ab initio simulationJ^. 

In the bulk, it is argued that liquid structure is de- 
termined to a significant degree by structures which are 
energetically locally favourabl^^. Such locally favoured 
structures can correspond to clusters which minimise the 
potential energy in isolation and may be catalogued and 
ordered by the number of particles they contaiiP^^'^. 
Structures which exhibit five-fold symmetry, in partic- 
ular, ar e very common and account for up to 2/3 of all 
particlej ^^*^ I As a result, the local dynamics can be 
determined to a significant degree by the geometry and 
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the symmetries of the locahy favoured structures. This 
observation is potentially significant to understand the 
glassy behaviour of supercooled liquids: if the symmetry 
of the locally favoured structure differs from that of the 
crystalline state, the path to crystallization is frustrated. 

It is thus natural to ask what happens to the con- 
centration of clusters near the liquid-gas interface. This 
will be a significant measure of local liquid structure. 
Two mechanisms will contribute to a significant change in 
cluster concentration near the interface. First, structures 
favoured by the bulk may be suppressed near the free sur- 
face. Second, clusters are non-isotropic, and might thus 
be ordered in a particular way with respect to the free 
surface. We will explore this ordering effect in the case 
of two very common clusters, following the nomencla- 
ture of Doye et. al.^^: the 5 A triangular bipyramid and 
7A pentagonal bipyramid, which are 3-and 5-membered 
rings with atoms attached above and below the plane of 
the ring. 

One might think that the free surface could lead to 
an enhancement of local fivefold symmetry. The reason 
is that in the case of deeply quenched liquids, the free 
interface can to an extent be thought of as a constrain- 
ing field. Now the structure induced in simple liquids by 
confinement such as hard walls has been shown, under 
the application of a second field, to induce some five-fold 
symmetr}^, which may be conjectured to relate to X-ray 
refiectometry experiments on lead which found evidence 
for local fivefold symmetry at the interface^^. Surpris- 
ingly, here our analysis suggests the opposite effect: 7 A 
pentagonal bipyramid clusters, which form the basic unit 
of five-fold symmetry in our methodolog}E^^, are sup- 
pressed disproportionately near the interface. 

To determine locally favoured structures, we use a 
novel method, the Topological Cluster Classification 
(TCC). This identifies small clusters of particles from 
within bulk phases which are topologically similar to a 
set of reference cluster In this case the reference struc- 
tures are formed by groups of 5 ^ m ^ 13 Morse par- 
ticles in isolation^^ and, depending on the range of the 
Morse potential, correspond to the global energy mini- 
mum clusters of the Lennard- Jones and Sodium poten- 
tials considered. We have recently explored the structure 
of the bulk Lennard- Jones liquid and a system similar to 
the Sodium model using this methocP"^. Contrary to the 
suggestion of Frank who conjectured that in the struc- 
ture of such liquids 13-membered icosahedral structures 
might be prevalent, we found that these only account 
for one particle in a thousand, according to our crite- 
ria. Instead five-fold symmetry stems from pentagonal 
bipyramids, which account for 54% of the particles at 
the Lennard- Jones triple point^^. 

This paper is organised as follows. First, we describe 
our methodology for identifying the interface, and local 
structure in section [ll| We then proceed to present and 



II. METHODOLOGY 



A. Models 



We use the Monte Carlo (MC) method to simulate a 
liquid in equilibrium with its vapour phase^^. We use 
the Lennard- Jones (LJ) 12-6 potential, where, for two 
particles separated by a distance r, the interaction energy 
U is given as 



(1) 



where f3 = l/ksT where A:^ is Boltzmann's constant and 
T is temperature. There are two parameters: slj is the 
strength of the attraction between the particles and cjlj 
determines the range of the interaction. From here on 
we quote all quantities in standard reduced units, where 
energies and lengths are normalized by £lj and cflj re- 
spectively. 

Now the triple point of Lennard- Jones lies around 
T^^"^ = 0.68, so the maximum degree of cooling rela- 
tive to the critical point is T^-/ /T^'^ = 0.63. We would 
like to realize even lower temperatures and this can be 
achieved by employing models with longer-ranged inter- 
actions which freeze at lower temperatures relative to 
criticality. In particular, the spherically symmetric model 
for Sodium (Na) introduced by Chacon et. al^ has a 
triple point around T^^ = 0.27 and critical point around 
j^Na _ ]^ 25, enabling a much larger degree of cooling, 
T^^^/T^^^ = 0.22. This model potential for Sodium reads 



(2) 



discuss our findings in section III before placing our find- 



ings in context in the concluding section IV 



where Aq = 437.96 eV, Ai = 0.18382 eV, a = 2.2322 
A"^ b = 0.2140 A"^ and Ri = 3.5344 A. A lengthscale 
o'Na = 3.48 A is defined as minimum of Ung along with a 
well depth SNa = UNai^'Na) = 0.1885 eV. Energies and 
lengths are then normalized by e^a and aNa respectively. 

In Fig. [T] we show the Lennard- Jones and Sodium po- 
tentials to be considered below. We also display the 
Morse potential, which will be used as a reference po- 
tential for the minimum energy clusters of the Sodium 
model. 

We truncate and shift the Lennard- Jones potential at 
finite range rcut = 2.5(T such that pULj{'2.5a) = 0. In the 
case of the model Sodium potential, we follow^ and trun- 
cate the model potential at Tcut = 2.5cr but, unlike the LJ 
case, we do not shift the potential, i.e. /3/7Ara(2.5<j) < 0. 



B. Simulation details 

We use Monte Carlo simulations in the NVT ensem- 
ble^^. The simulations consist of N = 16000 particles 
in cuboid boxes with sides LxI/x2L and with periodic 
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FIG. 1: Model potentials used. Blue line, Lennard- Jones 
(LJ); solid black line, Sodium (Na); dashed line, Morse po- 
tential with range parameter = 4.05. 



boundary conditions. This geometry ensures that, once 
equilibrated, the liquid forms a slab with two interfaces 
perpendicular to the z-direction. Assuming the interface 
does not overlap with itself, the z-position of the top or 
bottom interface may be expressed as a function of ix^y). 

We present results for two state points in the case of 
Lennard- Jones p = 0.32, T = 0.95 and p = 0.32, T = 
0.68, and in the case of Sodium p = 0.485, T = 0.542 
and p = 0.485, T = 0.270. These correspond to T/T^ = 
0.880, 0.630, 0.434 and 0.237 respectively where is the 
critical point of LJ or Na respectively, p is the overall 
density in the simulation box, which contains both liquid 
and vapour. This gives L = 29. 2a lj for the LJ state 
points and L = 25.baNa for the Sodium state points. 

For all simulations suitable liquid-vapour coexistence 
initial configurations are generated by first estimating ap- 
propriate bulk liquid pL and vapour pc densities at co- 
existence from the literature^ and then equilibrating 
a liquid slab of density Pl in coexistence with a vapour 
at pg for at least 300000 MC sweeps. Equilibration is 
ensured by both checking relaxation of the total poten- 
tial energy and by examining the evolution of the liquid- 
vapour density profile p{z). It has been shown that the 
coexistence densities and p{z) depend on the system size, 
the geometry of the box, and the boundary conditions of 
the simulatiorP^, as well as the truncation length and 
whether the potential is shifted^^. Once the density pro- 
file p{z) is seen to have stopped evolving, we say the 
system is at equilibrium and a production run of 10000 
MC sweeps is used to generate 100 independent config- 
urations for analysis. Four or more independent simula- 
tions are performed for each state point and the results 
averaged over all eight or more liquid-vapour interfaces 
obtained from these. 
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FIG. 2: Identification of the position of a capillary perturbed 
interface. The interface is shown in blue and the cube is 
split into columns The interface is assumed flat in 

each column and the position is extracted by a fitting step 
function (Eq. (|4|) to the columnar density profile pij{z). 
Note periodic boundaries apply in the x and y directions only. 

C. Determining the interface location 

To calculate the location of the interface we first fit 
p{z) with a hyperbolic tangent function in the following 
form 

.X PL^ pG ^ PL- pG , u f 
p{z) = + tanh i^---- j (3) 

where pL and pc are the bulk densities of the liquid and 
the gas at equilibrium coexistence, w quantifies the width 
of the interfacial region, and zq is the position of the inter- 
face with respect to the simulation coordinate axis z. As 
we will discuss below, w as determined by ([3| ignores the 
effect of capillary waves on the interface. Although it is 
difficult to accurately decouple the intrinsic width of the 
interface from that of a capillary broadened interfacd^, 
we show below, for all the state points studied here, ig- 
noring the effects of capillary broadening and using the 
fit ([3| to determine the position and orientation of the 
interface, that the conclusions drawn about the structure 
of the interface do not significantly change as a result of 
this approximation. 

To assess the effect of capillary waves, we note that 
Eq. ^ assumes that the interface position is indepen- 
dent of the local x, y position and is perfectly flat, which 
is clearly not the case if capillary waves are present. We 
imagine that the interfacial width may be decomposed 
into capillary-broadened and intrinsic components Wcw 
and Wi where w = Wcw + '^i and we are interested in the 
intrinsic width wj^. We determine the effect of surface 
roughening (capillary waves) as follows^^. Firstly we di- 
vide the simulation cube into columns of width and 



breadth ^ and height L for integer n as shown in Fig. 2 
Rather than Eq. ([3|, we model the density profile Pij(z) 
of column i,j with the step function 



\pG Z>Zo,ij, 



(4) 



where zo^ij is the position of the interface and i and j 
are integers denoting the column. We minimize the least 
square residuals between p^ji^z^ and Pstep 

(z) to find zo^ij. 

Re-scaling the z axis as pij{z — zo) gives columnar den- 
sity profiles centred on the position of the interface. Tak- 
ing the mean of these columnar profiles gives pav{z — zq) 
the averaged density profile with capillary waves at and 
above the columnar scale removed, which we then fit 
with the tanh function (|3| leading to a new interfacial 
width Wn- Due to capillary broadening, the square of the 
measured interfacial width has a logarithmic dependence 
on the column widtlP^, while extrapolating the column 
width to zero gives an estimate of the intrinsic width. 
However, for very small column widths, there are insuf- 
ficient statistics (coordinates) in the columns to reliably 
give the columnar density profiles Pij{z — zq). 

For the LJ T = 0.95 state point, where we expect the 
largest contribution Wcw from capillary waves, the mini- 
mum we obtain is Wn = 2.11{l)aLj^ compared to the fiat- 
interface value of w = 2.38(2)crLj. w as obtained from 
equation ([3| is only slightly (around 10%) larger than the 
minimum Wn obtained from the capillary wave analysis 
and therefore roughening of the interface is minimal over 
the lengthscale we measure. Following this analysis we 
conclude capillary waves make a rather limited impact 
on the interfacial widths we can measure and hereafter 
neglect the effect of capillary waves and simply fit the 
whole simulation box with Eq. ([3| to determine the lo- 
cation of the interface. We henceforth rescale our 2:-axis 
such that Zq = 0. 



D. Topological cluster classification (TCC) 

To analyse the structure, we first identify the bond 
network using a modified Voronoi construction with a 
maximum bond length Tcut = l-8cr and four-membered 
ring parameter fc = 0.82^^. Having identified the bond 
network, we use the topological cluster classification to 
determine the nature of the cluster. This analysis iden- 
tifies all the shortest path three, four and five membered 
rings in the bond network. These rings form the basic 
building blocks of our analysis. A 3-membered ring with 
atoms bound above and below is identified as a 5A tri- 
angular bipyramid, a 4-membered ring with two bound 
atoms as a 6A octahedron and a 5-membered ring as a 
7A pentagonal bipyramid, as shown in Fig. |3] The ad- 
ditional atoms bonded to the rings are termed spindle 
atoms. Now the ground state clusters for Lennard- Jones 
are known^^, but we are unaware of any work identify- 
ing ground state clusters for the Sodium model employed 
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FIG. 3: Structures of LJ and Na clusters, nomenclature fol- 
lowing^^, and their point group symmetries. The colours de- 
note the method used for cluster detection by the TCC algo- 
rithm: 3-, 4- and 5- membered ring atoms are grey, spindle 
atoms are yellow, additional atoms to a basic cluster are red 
(here just 8B which is detected from 7 A clusters), and rings 
are coloured white, blue, pink and green. 



here. However, ground state clusters for the variable 
range Morse potential have been calculated^. We there- 
fore assume that the clusters of the Morse potential with 
the relevant range (as discussed below) are also ground 
states for Sodium. We then use the TCC to find clusters 
which are global energy minima of the Lennard- Jones 
and Morse potentials, applying the latter to our Sodium 
results. Supporting this assumption is work which shows 
that for a model for Sodium with many-body interac- 
tions, the Gupta model, the ground state clusters are 
the same as those we infer for the Chacon and Tarazona 
model^^ 4^ 

The Morse potential is defined 



Po(o--r)j^gPo(cr-r) 



2), 



(5) 



where po is a range parameter and £m is the potential 
well depth. The Morse potential has a variable range and 
we choose the range such that it closely approximates 
the Sodium model. The extended law of corresponding 
states^^ provides a means by which different systems may 
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be compared with one another, by equating their reduced 
second virial coefficients 

B2 = (6) 

where aEFF is the effective hard sphere diameter and the 
second virial coefficient 

00 

^2 = 27r J drr'^ [1 - exp {-(3U{r))] . (7) 


The effective hard sphere diameter is defined as 
00 

(JEFF = J dr[l- exp {-f3UREp{r))] . (8) 


Varying the Morse range parameter such that 5| for 
Sodium and the Morse potential are equal gives po = 
4.05. Now the ground state minima for the Morse poten- 
tial are known^^ and here we assume those for po = 4.05 
also hold for UNa{^)- The potentials are compared in 
Fig. [1] We identify all topologically distinct Morse 
(po = 4.05) and Lennard- Jones clusters^^. In addition, 
we identify the FCC and HCP thirteen particle structures 
in terms of a central particle and its twelve nearest neigh- 
bours. We illustrate these clusters in Fig. [3] In the case 
of the Morse potential, for m > 7 there is more than one 
cluster which forms the ground state, depending on the 
range of the interaction^^. We therefore consider ground 
state clusters for each system and for m < 14. For more 
details sed^. 



III. RESULTS 

We begin our presentation of the results by showing 
the interfacial profiles of the different systems and state 
points. In Fig. [4] we see that the interfacial width drops 
with decreasing T/Tc. The lowest temperatures, and 
thus the smallest interfacial widths, are realized for the 
sodium system. We note the pronounced interfacial lay- 
ering at the lowest temperatur e!^ ' -'^^ l This indicates that 
the effect of the interface on the local ordering of the 
liquid extends far into the bulk in that case. 

We now proceed to the TCC analysis of the Lennard- 
Jones liquid at two different temperatures. We plot the 
total density of particles p{z) as well as the particle den- 
sities pc{z) for the most popular clusters. These are 
defined as the densities of all the particles belonging to 
a particular cluster. 

The densities are plotted on a logarithmic scale as a 
function of the distance z from the interface in Fig. |5j 
Note that the liquid side only is plotted > 0), and that, 
for Lennard- Jones, the density varies smoothly across the 
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FIG. 4: Interfacial profiles p(^z). From top, Lennard- Jones 
T = 0.95, T = 0.68, Sodium T = 0.542, T = 0.297. Data 
offset for clarity. The sodium interface at the lowest tempera- 
ture shows a very pronounced surface layering. The interfaces 
have all been centred at z — 0. 



interface for our treatmentP^. For each of the temper- 
atures, we also plot the cluster densities relative to the 
total density on a linear scale. The relative density is 
then defined as the density of particles in a given cluster 
pc{z) divided by the total density p{z). In the bulk it 
is seen that the most popular clusters, including the 7A 
pentagonal bipyramid, account for more than half of all 
the particles. 

As one moves towards the interface, the relative den- 
sity of clusters decreases. For the Lennard- Jones system, 
our analysis suggests that there is no enhancement of 
5-fold symmetry at this interface. In fact, as all clus- 
ter populations drop with the reducing density at the 
interface, we argue that there is more of a reduction in 
5-fold symmetry than anything else. In particular the 7 A 
pentagonal bipyramid cluster may be taken as a rough 
proxy for five-membered rings, the basic unit of 5-fold 
symmetry. These clearly follow the density profile, al- 
though the increase with depth into the liquid is much 
more marked than the overall density. Similar conclu- 
sions can be drawn about the 13A icosahedron. The de- 
crease in relative cluster density can be rationalized by 
observing that while potential energy considerations en- 
hance cluster populations^^, so too does packing^^. Thus 
we expect a decrease in cluster population at lower den- 
sity, for example upon approach to the interface. 

For example, for both T = 0.95 (T/T^^ = 0.880) and 
T^^/ = 0.68 {T/T^^ = 0.630), the relative decrease in 
density amounts to a factor of two in the case of 7A 
pentagonal bipyramid clusters. Since the change appears 
coupled to the overall density p, we see further change 
in the structure of the liquid near the surface for the 
temperatures studied. Note that although it would be 
attractive to investigate the effect of reducing the density 
on the cluster population in a bulk system of course this 
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FIG. 5: TCC analysis of the Lennard- Jones interface pc{z) denotes the density of particles belonging to a particular cluster 
type. p{z) reproduces the data in Fig. [i] (a) T = 0.95 {T/T^-^ = 0.880). (c) The triple point T^/ = 0.68 {T/T^-^ = 0.630). 
(d) Same data as (c) but normalised by p{z). Interface increases from 2; = 0, black interface density profile. Yellow shading is 
a guide to the eye, indicating the width of interfacial region. 



is not possible because the Lennard- Jones system would 
be unstable to vapour-liquid phase coexistence. 

To gain more insight into the effect of the interface 
on the cluster population, it would therefore be desirable 
to make the effect of the interface more pronounced by 
going to lower temperatures, thus making the interface 
sharper [for which we use the model Sodium potential 
Eq. (§]. 

Note that the triple point in general has a rather higher 
population of clusters than does T = 0.95, which is not 
unreasonable, as in addition to the higher density the en- 
hanced relative attractions at reduced temperature would 
be expected to promote clustering, which has been ob- 
served in the bulk^^. 

We next performed the same spatially resolved cluster 
analysis for Sodium at lower temperatures relative to crit- 
icality. We confirmed that, at higher T^^ = 0.63, the re- 
sults with Sodium were very similar to those for Lennard- 
Jones. As Fig. |4] has shown, the interface becomes nar- 
rower and at the lowest temperature there are oscillations 
of the bulk density profile indicating significant surface 
layering. This shows that the interface changes structur- 
ing into the liquid at longer ranges than in the case of 
Lennard- Jones. 

Once more we first show the absolute densities for all 
the clusters on a logarithmic scale, and then the densi- 



ties relative to the total on a linear scale (Fig. 6|. The 
latter are a direct measure of the influence of the inter- 
face on the local ordering. At the higher temperature, 
the decrease in the relative cluster density is more pro- 
nounced than before, but still remains modest. We note 
that once more there is certainly no increase of 7A pen- 
tagonal bipyramid cluster concentration owing to orien- 
tational ordering near the interface. 

At the lower temperature, 7A clusters account for 82% 
of the particles in the bulk, the relative density is re- 
duced to 31% near the interface, a somewhat larger drop 
in population. However, in stark contrast to the higher 
temperature state points, the effect on cluster concentra- 
tion extends far beyond the density decrease related to 
the interface, and into the bulk. This mirrors the layering 
as reflected in the oscillations of the density profile. The 
cluster concentrations are evidently a sensitive measure 
of the layering effect, and of the change in orient ational 
order it introduces. 

By contrast, the relative decrease of the triangular 
bipyramidal 5A clusters is relatively mild. As the clusters 
are smaller, this might be expected, since the perturba- 
tion by the interface on the compact 5A structure will 
be smaller. Other popular clusters, notably 6A octahe- 
dra and 8A D2d also show a similar behaviour. To gain 
more insight into the concentration of 5A and 7A clus- 
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FIG. 6: TCC analysis of the Sodium interface pc{z) denotes the density of particles belonging to a particular cluster type. 
p{z) reproduces the data in Fig. [i] (a) T = 0.542 (T/Tj^^ = 0.434). (c) Near the triple point T,^" = 0.297 (T/Tf " = 0.237). 
(b,d) Same data as (a,c) but normalised by p{z). Yellow shading is a guide to the eye and indicates the interfacial region. 



ters, we consider their orientational ordering near the 
interface. As seen in Figs. [3] and [7| both are viewed as a 
3-membered and 5-membered ring, respectively, with two 
spindle particles sticking out in a direction perpendicular 
(ignoring thermal fluctuations) to the ring. We neglect 
the 6A octahedron due to its high symmetry. To measure 
the orientation of these structures relative to the inter- 
face, we consider the angle between the spindle particles 
and the normal to the plane. From this we can construct 
an order parameter 1 /2 < 3 cos^ ^ — 1 > which is zero in 
the bulk, where all orientations are equal, unity in the 
case that all clusters are aligned with the axis perpendic- 
ular to the interface and —1/2 when clusters lie with the 
axis parallel to the interface. 

In Fig. [S] we plot the order parameter as function of 
z for both 5 A and 7A clusters for three different T/Tc 
ratios. Remarkably, the two types of clusters behave in 
an opposite way near the interface. While the 7 A clus- 
ters are aligned with the five-membered ring parallel to 
the interface, the 5 A clusters are aligned with its three 
membered ring perpendicular to the interface. This ef- 
fect increases strongly with decreasing temperature, so 
at that lowest temperature [Fig. (sjc)] in the case of 7A 
the angle is reduced, and the order parameter goes to 
1/2, half of what would be a perfectly ordered state. 

Thus it is likely that the 5-membered ring of the 7A 
clusters tends to be found in the plane of the interface. 
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FIG. 7: Definition of the orientational order parameter used 
in the case of the 5A triangular bipyramid and 7A pentagonal 
bipyramid. is the angle between the normal ruf to the plane 
of the interface (marked in blue) and the vector connecting 
the centres of the yellow spindle atoms. 



Assuming that 7A is in a configuration which minimises 
bond stretching/compression (i.e. zero-temperature po- 
tential energy minimum) , this would leave one end of the 
spindle sticking out from the free surface. This would 
be an energetically unfavourable situation, which could 
explain the suppression of 7A clusters near the interface. 



8 




0.4 

0.2 

0.0 h 
-0.2 
-0.4 




■7A 
■5A 



r/r "^'=0.434 



0.4 
0.2 
0.0 
-0.2 
-0.4 






2 


4 

z/o 


6 


1 


1 


' 1 


7 A - 

5A 


1 




1 


r/r ^'=0.237 " 

c 





2 


4 


6 



FIG. 8: The ordering of the 5 A and 7A clusters relative to 
interface, as function of the distance z from the interface. 
The angle is the angle between the vector that connects 
the spindle particles and the normal to the surface, (a) 
Lennard- Jones, T = 0.95, T/T^ = 0.880; (b) Sodium, T = 
0.542, T/T^ = 0.434; (c) Sodium, T = 0.297, T/T^ = 0.237. 



However the consequences of finite temperature may well 
play an important role here, a point to which we return 
in the next section. The situation is very different for 
5 A clusters, whose degree of ordering is also rather less 
pronounced. More crucially, though the ordering is such 
that the spindle remains in the fluid, while the three- 
membered ring is oriented parallel to the interface. This 
permits the cluster to remain inside the fluid relatively 



undistorted, and the effect on its concentration remains 
small. Fig. Sjc) also shows that the surface layering is 
well reflected by the orientational ordering of the clusters. 



IV. CONCLUSIONS 

We have shown that our topological cluster classifica- 
tion provides insight into the local liquid structure close 
to a free interface, and can directly probe for local struc- 
tures with fivefold symmetry. The analysis we have per- 
formed provides no evidence in support of an enhanced 
five-fold symmetry near the interface, however those five- 
fold symmetric structures that are present tend to be 
aligned with five-membered rings parallel to the inter- 
face. 

At higher temperatures, the predominant effect ap- 
pears to be related to the lowering in density induced 
by the interface. This is reflected by a drop in the pop- 
ulations of all the clusters we consider. Upon reducing 
temperature, the interfacial density profile exhibits lay- 
ering, and seems to induce a change in cluster population 
extending well into the dense liquid, long after the mean 
density is that of the bulk liquid. Thus we argue that our 
analysis indicates that five-fold symmetry is suppressed 
at the liquid- vapour interface in the systems we have con- 
sidered. 

We use our analysis to reveal orientational informa- 
tion on two basic clusters. Close to the interface, the 
five-fold symmetric 7A pentagonal bipyramid is oriented 
with its five-membered ring lying parallel to the interface, 
conversely the 5A triangular bipyramid is oriented with 
its 3-membered ring perpendicular to the interface. Both 
these effects are strongly enhanced at lower temperature. 

Our analysis is topological in nature. That is to say, 
the structures identified are based on bonds, and we do 
not consider structural distortions. In particular, we note 
that the 7A pentagonal bipyramid, aligned with a sur- 
face, might leave a spindle atom exposed. However, for 
a small amount of strain, this exposed spindle could in 
fact lie very close to the plane of the five-membered ring. 
We suspect that it may do so at finite temperature. 

We identify two possible extensions of this work. It 
would be interesting to consider the effect of a wall or 
similar external field, as an alternative to the free in- 
terface. It would also be interesting to investigate the 
behaviour of clusters other than those which are min- 
imum energy ground states in isolation, for example 5- 
membered rings. We have seen that the population of 7 A 
pentagonal bipyramids is reduced near the interface, the 
topological basis of our analysis precludes ruling out that 
an exposed spindle atom might be related to this suppres- 
sion, although at finite temperature we expect this effect 
to be small. An extension of our analysis to consider dis- 
tortions in the clusters identified would be helpful. Both 
these possibilities will be investigated in the future. 
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